US COVID Data

Data Prep

us <- read.csv("https://raw.githubusercontent.com/JaclynCoate/6373_Time_Series/master/TermProject/Data/USdaily.7.26.20.csv", header = T, strip.white = T)
us <- transform(us, date = as.Date(as.character(date), "%Y%m%d"))
us <- subset(us, select = -c(states, dateChecked, hospitalized, lastModified, total, posNeg, totalTestResultsIncrease, hash))
us[is.na(us)] <- 0
us = us[order(as.Date(us$date, format = "%Y%m%d")),]
head(us)
##           date positive negative pending hospitalizedCurrently
## 187 2020-01-22        2        0       0                     0
## 186 2020-01-23        2        0       0                     0
## 185 2020-01-24        2        0       0                     0
## 184 2020-01-25        2        0       0                     0
## 183 2020-01-26        2        0       0                     0
## 182 2020-01-27        2        0       0                     0
##     hospitalizedCumulative inIcuCurrently inIcuCumulative
## 187                      0              0               0
## 186                      0              0               0
## 185                      0              0               0
## 184                      0              0               0
## 183                      0              0               0
## 182                      0              0               0
##     onVentilatorCurrently onVentilatorCumulative recovered death
## 187                     0                      0         0     0
## 186                     0                      0         0     0
## 185                     0                      0         0     0
## 184                     0                      0         0     0
## 183                     0                      0         0     0
## 182                     0                      0         0     0
##     totalTestResults deathIncrease hospitalizedIncrease negativeIncrease
## 187                2             0                    0                0
## 186                2             0                    0                0
## 185                2             0                    0                0
## 184                2             0                    0                0
## 183                2             0                    0                0
## 182                2             0                    0                0
##     positiveIncrease
## 187                0
## 186                0
## 185                0
## 184                0
## 183                0
## 182                0

Stationarity: Current Hospitalizations

It is difficult to assume stationarity for this data due to multiple factors. We are working under the assumption that COVID is a novel virus and cases as well as hospitalizations will eventually return to zero. This being said our current modeling techniques do things such as return to the median or mimic the previously seen trends. Also, we see a heavy wandering trend in both new cases and hospitalization would be dependent on this as well as time. We will review the data and see what, if any, non-stationary components reveal themselves and model the data accordingly.

Current Hospitalization Realization

Traits:

  • Heavy wandering behavior
  • What appears to be some noise that could be pseudo-cyclic behavior hidden by the large numbers.
ggplot(data = us, aes(x=date, y=hospitalizedCurrently))+
  geom_line(color="orange")+
  labs(title = "Current COVID Hospitalized Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

Independent Variable Review

When reviewing the variables and the scatter plot from our previous EDA we can see that there are some correlations that were expected, for example currentHospitalizations is correlated with the variables that refelct ICU and Ventilaor patients. These metrics wouldn’t exist without hospitalizations. These variables also cannot help up predict hospitalizations because they after occurances. So we will actually be leveraing variables such as positive increase in order to see if there is some sort of correlation between hospitalized patients and the number of positive cases.

  • Positive Increase Trend
ggplot(data = us, aes(x=date, y=positiveIncrease))+
  geom_line(color="orange")+
  labs(title = "Increase in COVID Cases US", y = "Thousands", x = "") +
    theme_fivethirtyeight()

Multivariate MLR w/ Correlated Errors for Currently Hospitalized Patients

Forecast Independent Variables: Increase Positive Cases

  1. Evaluation: Slowly dampening ACF and heavy wandering consistent with a (1-B) factor. As well as frequency peaks at 0 and .14. Consistent with (1-B) and seasonailty component of 7.
#a
plotts.sample.wge(us$positiveIncrease)

## $autplt
##  [1] 1.0000000 0.9740745 0.9439116 0.9109437 0.8869924 0.8711719 0.8605261
##  [8] 0.8462935 0.8144161 0.7772852 0.7356766 0.7045767 0.6829609 0.6646632
## [15] 0.6432819 0.6088202 0.5688872 0.5293417 0.5001353 0.4761589 0.4592724
## [22] 0.4426271 0.4148783 0.3789333 0.3458221 0.3182895
## 
## $freq
##  [1] 0.005347594 0.010695187 0.016042781 0.021390374 0.026737968
##  [6] 0.032085561 0.037433155 0.042780749 0.048128342 0.053475936
## [11] 0.058823529 0.064171123 0.069518717 0.074866310 0.080213904
## [16] 0.085561497 0.090909091 0.096256684 0.101604278 0.106951872
## [21] 0.112299465 0.117647059 0.122994652 0.128342246 0.133689840
## [26] 0.139037433 0.144385027 0.149732620 0.155080214 0.160427807
## [31] 0.165775401 0.171122995 0.176470588 0.181818182 0.187165775
## [36] 0.192513369 0.197860963 0.203208556 0.208556150 0.213903743
## [41] 0.219251337 0.224598930 0.229946524 0.235294118 0.240641711
## [46] 0.245989305 0.251336898 0.256684492 0.262032086 0.267379679
## [51] 0.272727273 0.278074866 0.283422460 0.288770053 0.294117647
## [56] 0.299465241 0.304812834 0.310160428 0.315508021 0.320855615
## [61] 0.326203209 0.331550802 0.336898396 0.342245989 0.347593583
## [66] 0.352941176 0.358288770 0.363636364 0.368983957 0.374331551
## [71] 0.379679144 0.385026738 0.390374332 0.395721925 0.401069519
## [76] 0.406417112 0.411764706 0.417112299 0.422459893 0.427807487
## [81] 0.433155080 0.438502674 0.443850267 0.449197861 0.454545455
## [86] 0.459893048 0.465240642 0.470588235 0.475935829 0.481283422
## [91] 0.486631016 0.491978610 0.497326203
## 
## $db
##  [1]  14.3997611  15.8955909   7.7642205   8.7050186   3.8198773
##  [6]  -1.1863576   2.6815084  -0.4667544  -1.0582274  -3.9729372
## [11]  -4.4388886  -5.3739230  -7.0574666  -3.9001118  -5.3732547
## [16]  -5.9203161  -6.0265835  -6.3073466  -8.7427682  -8.5315257
## [21]  -8.1694048  -8.8726978  -6.2228751  -5.6255282  -5.7869081
## [26]  -2.3190255  -0.6895940 -18.2510780 -12.7252753 -17.9895004
## [31] -13.1495229 -17.1809062 -12.9879824 -18.5885318 -14.0950716
## [36] -19.2823521 -16.1802599 -22.2193540 -17.1778910 -15.9864725
## [41] -13.8720977 -14.3606703 -12.3482825 -19.0159734 -12.8299009
## [46] -16.1634117 -19.2402721 -16.9940391 -19.9193164 -24.4462952
## [51] -12.2184745 -14.1347405 -11.1407485 -19.5608640 -22.3406178
## [56] -26.2948840 -15.9219942 -19.7896219 -16.3429699 -16.6545137
## [61] -13.2991307 -20.9465604 -12.5392646 -23.3015328 -21.5282804
## [66] -17.4748203 -26.6390880 -24.5923917 -17.5647166 -19.7282494
## [71] -17.6845420 -23.0825883 -19.0743556 -19.4958696 -19.5138333
## [76] -20.4270518 -18.7343253 -16.8982195 -14.1306252 -14.2054286
## [81] -15.0466299 -13.4474119 -14.7662082 -15.2848488 -14.7767269
## [86] -17.2543647 -17.5524024 -16.7680943 -23.6931393 -19.6888637
## [91] -22.9328710 -17.1805763 -15.7422487
## 
## $dbz
##  [1]  12.1744244  11.8068871  11.1912920  10.3234897   9.1987290
##  [6]   7.8134517   6.1690554   4.2794940   2.1855400  -0.0227779
## [11]  -2.1856123  -4.0930116  -5.5814407  -6.6274820  -7.3078903
## [16]  -7.6988921  -7.8526856  -7.8271468  -7.6919809  -7.5061868
## [21]  -7.3030888  -7.0977786  -6.9055053  -6.7545171  -6.6859448
## [26]  -6.7448940  -6.9710578  -7.3933729  -8.0283548  -8.8797888
## [31]  -9.9376266 -11.1745795 -12.5395304 -13.9486256 -15.2804168
## [36] -16.3913300 -17.1665125 -17.5812779 -17.7097183 -17.6695854
## [41] -17.5637030 -17.4562894 -17.3753921 -17.3239844 -17.2913960
## [46] -17.2628022 -17.2262162 -17.1767997 -17.1186074 -17.0640034
## [51] -17.0309217 -17.0382729 -17.1002705 -17.2209453 -17.3905128
## [56] -17.5857017 -17.7761429 -17.9367144 -18.0606665 -18.1649087
## [61] -18.2832201 -18.4521578 -18.6979308 -19.0286543 -19.4319317
## [66] -19.8765383 -20.3176913 -20.7053075 -20.9929698 -21.1442791
## [71] -21.1356798 -20.9585908 -20.6228760 -20.1586258 -19.6118560
## [76] -19.0349178 -18.4769322 -17.9782467 -17.5690679 -17.2703544
## [81] -17.0951532 -17.0493949 -17.1318377 -17.3332302 -17.6350503
## [86] -18.0085946 -18.4157858 -18.8133902 -19.1613062 -19.4324915
## [91] -19.6189508 -19.7292633 -19.7788462
  1. Differencing Data: much more stationary data and have surfaced a seasonaly component = 7 seen on ACF peaks at 7, 14, 21
#b
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

  1. Seaonality Transformation: stationary data, and an ACF that reflects data other than whitenoise to be modeled.
#c
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

  1. Model IDing: diagnose with aic.wge to determine phi, thetas: Both AIC/BIC select ARMA(4,2)
#d
aic5.wge(inpos.diff.seas)
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  aic
##       p    q        aic
## 15    4    2   15.62368
## 9     2    2   15.74701
## 18    5    2   15.74771
## 12    3    2   15.75351
## 7     2    0   15.75551
aic5.wge(inpos.diff.seas, type = "bic")
## ---------WORKING... PLEASE WAIT... 
## 
## 
## Five Smallest Values of  bic
##       p    q        bic
## 15    4    2   15.74832
## 2     0    1   15.79201
## 7     2    0   15.80893
## 4     1    0   15.81457
## 3     0    2   15.81522
  1. White Noise Test: Reject the H_null and accept this data is not white noise
#e
acf(inpos.diff.seas)

ljung.wge(inpos.diff.seas)$pval
## Obs -0.3711871 -0.02127283 0.07854862 0.05194601 -0.1044838 0.3276797 -0.4658847 0.1565306 0.1284186 -0.1060791 -0.07693327 0.1030723 -0.06391648 0.05002847 0.004471288 -0.03640442 -0.01551996 0.04133725 -0.01110058 -0.08153895 0.06570572 -0.03709739 -0.02390062 0.04588417
## [1] 1.246447e-12
  1. Estiamte Phis, Thetas
  2. Model Building
#f
est.inpos.seas = est.arma.wge(inpos.diff.seas, p = 4, q = 2)
## 
## Coefficients of Original polynomial:  
## -1.8727 -1.4974 -0.4735 0.0283 
## 
## Factor                 Roots                Abs Recip    System Freq 
## 1+1.1293B+0.6983B^2   -0.8086+-0.8822i      0.8357       0.3681
## 1+0.7944B             -1.2588               0.7944       0.5000
## 1-0.0510B              19.6244               0.0510       0.0000
##   
## 
mean(us$positiveIncrease)
## [1] 22567.12
#g
  1. Forecast
  • 7 Day
#7 day
inpos.preds7 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 7)

  • 90 Day
#90 day
inpos.preds90 = fore.aruma.wge(us$positiveIncrease, phi = est.inpos.seas$phi, theta = est.inpos.seas$theta, d=1, s=7, n.ahead = 90)

  1. Plotting forecasts
  • 7 Day Forecast
#7 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,195), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "7 Day Increase in COVID Cases Forecast")
lines(seq(188, 194,1), inpos.preds7$f, type = "l", col = "red")

  • 90 Day Forecast
#90 day forecast
plot(seq(1,187,1), us$positiveIncrease, type = "l", xlim = c(0,280), ylim = c(0,80000), ylab = "Increase in COVID Cases", main = "90 Day Increase in COVID Cases Forecast")
lines(seq(188, 277,1), inpos.preds90$f, type = "l", col = "red")

Forecast Dependent Variable: MLR w/ Correlated Errors for Currently Hospitalized Patients using Positive Increase

  1. Data prep and recognizing differencing and seaonal components of Currently Hospitalized
#Selecting only those dates with reported current hospitilizations
us <- dplyr::slice(us,56:187)
invisible(us)
#Stationarize Currently Hospitalized Variable
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

acf(us.diff)
us.diff.seas = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

acf(us.diff.seas)
#Stationarize Increase Positive Variable
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

acf(inpos.diff)
inpos.diff.seas = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

acf(inpos.diff.seas)

  1. Fit simple regression model predicting hospitalizedCurrently with positiveIncrease.
mlr.fit = lm(us.diff.seas~inpos.diff.seas, data=cbind(data.frame(us.diff.seas), data.frame(inpos.diff.seas)))
summary(mlr.fit)
## 
## Call:
## lm(formula = us.diff.seas ~ inpos.diff.seas, data = cbind(data.frame(us.diff.seas), 
##     data.frame(inpos.diff.seas)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7149.3  -554.2    73.4   611.2  6506.7 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)     -13.32614  143.23299  -0.093  0.92603   
## inpos.diff.seas   0.11733    0.04227   2.776  0.00637 **
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1595 on 122 degrees of freedom
## Multiple R-squared:  0.0594, Adjusted R-squared:  0.0517 
## F-statistic: 7.705 on 1 and 122 DF,  p-value: 0.006375
AIC(mlr.fit)
## [1] 2184.712
acf(mlr.fit$residuals)

plot(mlr.fit$residuals)

  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
  • Below we can see that our aic.wge function has selected an ARMA(1,2) model for modeling our cmort information.
mlr.phis= aic.wge(mlr.fit$residuals)
mlr.phis
## $type
## [1] "aic"
## 
## $value
## [1] 14.53403
## 
## $p
## [1] 5
## 
## $q
## [1] 2
## 
## $phi
## [1] -0.2465069 -0.4219100  0.3119255  0.2116968  0.2024803
## 
## $theta
## [1] -0.4657537 -0.9566014
## 
## $vara
## [1] 1803071
  1. Now forcast with arima function with phi’s from above coefficients and ARMA(1,2) model,
mlr.fit.arima = arima(us$hospitalizedCurrently, order = c(5,1,2), seasonal = list(order = c(1,0,0), period = 7), xreg = us$positiveIncrease)
## Warning in log(s2): NaNs produced
AIC(mlr.fit.arima)
## [1] 2224.09
  1. Run test to see if data is white noise: confirmed white noise, continue forward
acf(mlr.fit.arima$resid) 

ltest1 = ljung.wge(mlr.fit.arima$resid) 
## Obs -0.006649913 -0.009334226 -0.01070918 0.01827832 -0.02771128 -0.01061308 -0.01123975 0.08781402 -0.05809279 0.07500195 0.05661161 -0.08527475 -0.007553036 0.1098068 0.002427647 -0.010962 -0.04242877 -0.001526285 -0.05612271 -0.03313627 0.02930685 -0.007307643 -0.06056842 -0.03342384
ltest1$pval
## [1] 0.999208
ltest2 = ljung.wge(mlr.fit.arima$resid, K= 48)
## Obs -0.006649913 -0.009334226 -0.01070918 0.01827832 -0.02771128 -0.01061308 -0.01123975 0.08781402 -0.05809279 0.07500195 0.05661161 -0.08527475 -0.007553036 0.1098068 0.002427647 -0.010962 -0.04242877 -0.001526285 -0.05612271 -0.03313627 0.02930685 -0.007307643 -0.06056842 -0.03342384 0.08407797 -0.006875273 -0.002312484 0.06777622 -0.09591675 -0.06198868 -0.01331578 -0.03441461 -0.04091228 0.0270524 0.06615092 0.06538511 -0.09261011 -0.04939741 -0.03187124 -0.04934637 -0.01834029 -0.02031817 -0.001873916 0.02337696 -0.004828854 0.006148926 -0.002455308 0.004810674
ltest2$pval
## [1] 0.9999866
  1. Load the forecasted increase positive in a data frame
  • 7 Day
#7 Day Case Increase
regs7 = data.frame(positiveIncrease = inpos.preds7$f)
invisible(regs7)
  • 90 Day
#90 Day Case Increase
regs90 = data.frame(positiveIncrease = inpos.preds90$f)
invisible(regs90)
  1. Predictions
  • 7 Day
mlr1.preds7 = predict(mlr.fit.arima, newxreg = regs7, n.ahead = 7)
invisible(mlr1.preds7)
  • 90 Day
mlr1.preds90 = predict(mlr.fit.arima, newxreg = regs90, n.ahead =90)
invisible(mlr1.preds90)
  1. Plotted Forecasts
  • 7 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,140), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,139,1), mlr1.preds7$pred, type = "l", col = "red")

  • 90 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,223), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "90 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,222,1), mlr1.preds90$pred, type = "l", col = "red")

  1. Windowed ASE

Forecast Dependent Variable: MLR w/ Correlated Errors for Currently Hospitalized Patients w/ Positive Increase & Trend

  1. Fit simple regression model predicting hospitalizedCurrently with positiveIncrease and trend
#creating trend
time <- seq(1,124,1)
#fitting model
mlr.fit.t = lm(us.diff.seas~inpos.diff.seas+time, data=cbind(data.frame(us.diff.seas), data.frame(inpos.diff.seas)))
summary(mlr.fit.t)
## 
## Call:
## lm(formula = us.diff.seas ~ inpos.diff.seas + time, data = cbind(data.frame(us.diff.seas), 
##     data.frame(inpos.diff.seas)))
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -7111.0  -524.5    73.9   617.7  6540.6 
## 
## Coefficients:
##                  Estimate Std. Error t value Pr(>|t|)   
## (Intercept)      32.35112  289.27171   0.112  0.91114   
## inpos.diff.seas   0.11724    0.04244   2.762  0.00663 **
## time             -0.73096    4.01658  -0.182  0.85590   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1601 on 121 degrees of freedom
## Multiple R-squared:  0.05966,    Adjusted R-squared:  0.04412 
## F-statistic: 3.839 on 2 and 121 DF,  p-value: 0.02419
AIC(mlr.fit.t)
## [1] 2186.678
acf(mlr.fit.t$residuals)

plot(mlr.fit.t$residuals)

  1. Diagnose with aic.wge and see what the p, q would be with residuals from above model
    1. Below we can see that our aic.wge function has selected an ARMA(5,2) model for modeling our cmort information.
mlr.phis= aic.wge(mlr.fit.t$residuals)
mlr.phis
## $type
## [1] "aic"
## 
## $value
## [1] 14.53329
## 
## $p
## [1] 5
## 
## $q
## [1] 2
## 
## $phi
## [1] -0.2466126 -0.4220817  0.3118538  0.2119376  0.2027949
## 
## $theta
## [1] -0.4655288 -0.9564592
## 
## $vara
## [1] 1801724
  1. Now forcast with arima function with phi’s from above coefficients and ARMA(1,2) model,
Time <- seq(1,132,1)
mlr.fit.arima.t = arima(us$hospitalizedCurrently, order = c(5,1,2), seasonal = list(order = c(1,0,0), period = 7), xreg = cbind(us$positiveIncrease, Time))
AIC(mlr.fit.arima.t)
## [1] 2230.756
  1. Run test to see if data is white noise; confirmed white noise continue forward.
acf(mlr.fit.arima.t$resid) 

ltest1 = ljung.wge(mlr.fit.arima.t$resid) 
## Obs 0.002267493 0.004743329 0.004085712 0.008584111 0.003613306 0.03592544 -0.02552583 0.08439651 -0.1060931 0.01276058 0.02640347 -0.06061388 0.05268146 0.1458814 0.008052735 -0.04255608 -0.0935536 -0.04208266 -0.0498013 -0.001717023 0.05934226 0.01341049 -0.07137289 -0.07279666
ltest1$pval
## [1] 0.9823161
ltest2 = ljung.wge(mlr.fit.arima.t$resid, K= 48)
## Obs 0.002267493 0.004743329 0.004085712 0.008584111 0.003613306 0.03592544 -0.02552583 0.08439651 -0.1060931 0.01276058 0.02640347 -0.06061388 0.05268146 0.1458814 0.008052735 -0.04255608 -0.0935536 -0.04208266 -0.0498013 -0.001717023 0.05934226 0.01341049 -0.07137289 -0.07279666 0.05320601 -0.01313386 0.01794322 0.1027134 -0.07595104 -0.06199621 -0.03306343 -0.06206906 -0.03944973 0.03988941 0.07572305 0.07297186 -0.1029145 -0.06857141 -0.04403521 -0.04289784 0.005428483 -0.003871671 -7.833593e-05 0.005208028 -0.03085214 -0.01036108 0.00645461 0.02868683
ltest2$pval
## [1] 0.9990337
  1. Load the forecasted increase positive in a data frame
    1. 7 Day
#7 Day Case Increase
regs7t = data.frame(cbind(positiveIncrease = inpos.preds7$f, Time = seq(133,139)))
invisible(regs7)
b. 90 Day
#90 Day Case Increase
regs90t = data.frame(cbind(positiveIncrease = inpos.preds90$f, Time = seq(133,222)))
invisible(regs90)
  1. Predictions
    1. 7 Day
mlr1.preds7.t = predict(mlr.fit.arima.t, newxreg = regs7t, n.ahead = 7)
invisible(mlr1.preds7.t)
b. 90 Day
mlr1.preds90.t = predict(mlr.fit.arima.t, newxreg = regs90t, n.ahead =90)
invisible(mlr1.preds90.t)
  1. Plotted Forecasts
    1. 7 Day
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,140), ylim = c(0,60000), ylab = "Currently Hospitalized COVID Cases", main = "7 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,139,1), mlr1.preds7.t$pred, type = "l", col = "red")

b. 90 Day

plot(seq(1,132,1), us$hospitalizedCurrently, type = "l",  xlim = c(0,223), ylim = c(0,85000), ylab = "Currently Hospitalized COVID Cases", main = "90 Day Forecast for COVID Hospitalized Cases")
lines(seq(133,222,1), mlr1.preds90.t$pred, type = "l", col = "red")

  1. Windowed ASE
  1. 7 Day

  2. 90 Day

Forecast Dependent Variable: MLR w/ Correlated Errors (lagged variables) for Currently Hospitalized Patients w/ Positive Increase & Trend

With a quick check, we can see that there is no lag correlationed between the increase of COVID patients and hospitalized patients, like we theorized there might be. So we will not model an MLR w/ correlated errors on lagged variables.

ccf(us$positiveIncrease,us$hospitalizedCurrently)

Multivariate VAR Model

Since we are working with a seasonal and transformation component but it requires an additional transformation for the total positive COVID cases to become stationary for evaluation, we will only use positive increase of cases to predict currently hospitzliaed cases for the VAR model.

  1. Differenced and Seasonal Transformation VAR Model
#Positive Cases Transformations
#pos.diff = artrans.wge(us$positive, phi.tr = 1, lag.max = 100)
#pos.diff.2 = artrans.wge(pos.diff, phi.tr = 1, lag.max = 100)
#pos.trans = artrans.wge(pos.diff.2,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)
#Increase Positive Transformation
inpos.diff = artrans.wge(us$positiveIncrease, phi.tr = 1, lag.max = 100)

inpos.trans = artrans.wge(inpos.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

#Current Hospitalization Transformations
us.diff = artrans.wge(us$hospitalizedCurrently, phi.tr = 1, lag.max = 100)

currhosp.trans = artrans.wge(us.diff,phi.tr = c(0,0,0,0,0,0,1), lag.max = 100)

  1. Use Var select to find best model and fit model: p = 8 for lowest AIC
#VARselect to choose lag
VARselect(cbind(currhosp.trans, inpos.trans), lag.max = 10, type= "both")
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      7      7      2      7 
## 
## $criteria
##                   1            2            3            4            5
## AIC(n) 3.093857e+01 3.083329e+01 3.079091e+01 3.076976e+01 3.069554e+01
## HQ(n)  3.101650e+01 3.095018e+01 3.094676e+01 3.096458e+01 3.092932e+01
## SC(n)  3.113058e+01 3.112131e+01 3.117493e+01 3.124979e+01 3.127158e+01
## FPE(n) 2.731960e+13 2.459304e+13 2.357880e+13 2.309552e+13 2.145770e+13
##                   6            7            8            9           10
## AIC(n) 3.066656e+01 3.051427e+01 3.053424e+01 3.059674e+01 3.059649e+01
## HQ(n)  3.093930e+01 3.082598e+01 3.088492e+01 3.098638e+01 3.102509e+01
## SC(n)  3.133860e+01 3.128233e+01 3.139830e+01 3.155681e+01 3.165257e+01
## FPE(n) 2.086403e+13 1.793906e+13 1.833019e+13 1.955149e+13 1.959499e+13
#fit model
usdiff.var.fit = VAR(cbind(currhosp.trans, inpos.trans), type = "both", p = 2)
  1. Predictions for Difference
    1. 7 Day
#7 day predictions
pred.var7 = predict(usdiff.var.fit, n.ahead = 7)
pred.var7$fcst$currhosp.trans
##            fcst     lower    upper       CI
## [1,]  -42.12350 -2942.410 2858.163 2900.286
## [2,] -385.17371 -3346.493 2576.145 2961.319
## [3,]  -80.65766 -3262.529 3101.214 3181.871
## [4,] -124.14161 -3327.324 3079.041 3203.182
## [5,]  -46.62794 -3284.538 3191.282 3237.910
## [6,]  -43.58380 -3288.207 3201.039 3244.623
## [7,]  -11.23201 -3262.277 3239.813 3251.045
b. 90 Day
#90 day predictions
pred.var90 = predict(usdiff.var.fit, n.ahead = 90)
pred.var90$fcst$currhosp.trans
##              fcst     lower    upper       CI
##  [1,]  -42.123501 -2942.410 2858.163 2900.286
##  [2,] -385.173707 -3346.493 2576.145 2961.319
##  [3,]  -80.657665 -3262.529 3101.214 3181.871
##  [4,] -124.141611 -3327.324 3079.041 3203.182
##  [5,]  -46.627942 -3284.538 3191.282 3237.910
##  [6,]  -43.583795 -3288.207 3201.039 3244.623
##  [7,]  -11.232008 -3262.277 3239.813 3251.045
##  [8,]   -6.349200 -3259.226 3246.528 3252.877
##  [9,]    7.305513 -3246.870 3261.481 3254.175
## [10,]   12.421029 -3242.210 3267.053 3254.632
## [11,]   18.773716 -3236.131 3273.679 3254.905
## [12,]   22.480950 -3232.533 3277.495 3255.014
## [13,]   26.203504 -3228.870 3281.277 3255.073
## [14,]   28.923414 -3226.175 3284.022 3255.099
## [15,]   31.503265 -3223.608 3286.615 3255.112
## [16,]   33.699320 -3221.418 3288.817 3255.117
## [17,]   35.776726 -3219.344 3290.897 3255.120
## [18,]   37.697587 -3217.424 3292.819 3255.122
## [19,]   39.549962 -3215.572 3294.672 3255.122
## [20,]   41.334045 -3213.789 3296.457 3255.123
## [21,]   43.082251 -3212.041 3298.205 3255.123
## [22,]   44.799892 -3210.323 3299.923 3255.123
## [23,]   46.499486 -3208.623 3301.622 3255.123
## [24,]   48.185069 -3206.938 3303.308 3255.123
## [25,]   49.861824 -3205.261 3304.985 3255.123
## [26,]   51.532055 -3203.591 3306.655 3255.123
## [27,]   53.198022 -3201.925 3308.321 3255.123
## [28,]   54.860928 -3200.262 3309.984 3255.123
## [29,]   56.521787 -3198.601 3311.645 3255.123
## [30,]   58.181203 -3196.942 3313.304 3255.123
## [31,]   59.839642 -3195.283 3314.963 3255.123
## [32,]   61.497398 -3193.626 3316.620 3255.123
## [33,]   63.154688 -3191.968 3318.278 3255.123
## [34,]   64.811654 -3190.311 3319.935 3255.123
## [35,]   66.468399 -3188.655 3321.591 3255.123
## [36,]   68.124990 -3186.998 3323.248 3255.123
## [37,]   69.781476 -3185.341 3324.904 3255.123
## [38,]   71.437889 -3183.685 3326.561 3255.123
## [39,]   73.094252 -3182.029 3328.217 3255.123
## [40,]   74.750580 -3180.372 3329.874 3255.123
## [41,]   76.406884 -3178.716 3331.530 3255.123
## [42,]   78.063172 -3177.060 3333.186 3255.123
## [43,]   79.719449 -3175.404 3334.842 3255.123
## [44,]   81.375717 -3173.747 3336.499 3255.123
## [45,]   83.031981 -3172.091 3338.155 3255.123
## [46,]   84.688240 -3170.435 3339.811 3255.123
## [47,]   86.344498 -3168.778 3341.467 3255.123
## [48,]   88.000753 -3167.122 3343.124 3255.123
## [49,]   89.657007 -3165.466 3344.780 3255.123
## [50,]   91.313260 -3163.810 3346.436 3255.123
## [51,]   92.969513 -3162.153 3348.092 3255.123
## [52,]   94.625766 -3160.497 3349.749 3255.123
## [53,]   96.282018 -3158.841 3351.405 3255.123
## [54,]   97.938270 -3157.185 3353.061 3255.123
## [55,]   99.594521 -3155.528 3354.717 3255.123
## [56,]  101.250773 -3153.872 3356.374 3255.123
## [57,]  102.907025 -3152.216 3358.030 3255.123
## [58,]  104.563276 -3150.560 3359.686 3255.123
## [59,]  106.219528 -3148.903 3361.342 3255.123
## [60,]  107.875779 -3147.247 3362.999 3255.123
## [61,]  109.532031 -3145.591 3364.655 3255.123
## [62,]  111.188282 -3143.935 3366.311 3255.123
## [63,]  112.844534 -3142.278 3367.967 3255.123
## [64,]  114.500785 -3140.622 3369.624 3255.123
## [65,]  116.157037 -3138.966 3371.280 3255.123
## [66,]  117.813288 -3137.310 3372.936 3255.123
## [67,]  119.469540 -3135.653 3374.592 3255.123
## [68,]  121.125791 -3133.997 3376.249 3255.123
## [69,]  122.782043 -3132.341 3377.905 3255.123
## [70,]  124.438294 -3130.685 3379.561 3255.123
## [71,]  126.094546 -3129.028 3381.218 3255.123
## [72,]  127.750797 -3127.372 3382.874 3255.123
## [73,]  129.407049 -3125.716 3384.530 3255.123
## [74,]  131.063300 -3124.060 3386.186 3255.123
## [75,]  132.719552 -3122.403 3387.843 3255.123
## [76,]  134.375803 -3120.747 3389.499 3255.123
## [77,]  136.032055 -3119.091 3391.155 3255.123
## [78,]  137.688306 -3117.435 3392.811 3255.123
## [79,]  139.344558 -3115.778 3394.468 3255.123
## [80,]  141.000809 -3114.122 3396.124 3255.123
## [81,]  142.657061 -3112.466 3397.780 3255.123
## [82,]  144.313312 -3110.810 3399.436 3255.123
## [83,]  145.969564 -3109.153 3401.093 3255.123
## [84,]  147.625815 -3107.497 3402.749 3255.123
## [85,]  149.282067 -3105.841 3404.405 3255.123
## [86,]  150.938318 -3104.185 3406.061 3255.123
## [87,]  152.594570 -3102.528 3407.718 3255.123
## [88,]  154.250821 -3100.872 3409.374 3255.123
## [89,]  155.907073 -3099.216 3411.030 3255.123
## [90,]  157.563324 -3097.560 3412.686 3255.123
  1. Calculate Actual Forecasts from Predicted Differences
    1. 7 Day
startingPoints7 = us$hospitalizedCurrently[126:132]
currHospForecasts7 = pred.var7$fcst$currhosp.trans[,1:3] + startingPoints7
b. 90 Day
startingPoints90 = us$hospitalizedCurrently[43:132]
currHospForecasts90 = pred.var90$fcst$currhosp.trans[,1:3] + startingPoints90
  1. Plotting Forecasts
    1. 7 Day
#7 day Forecasts
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l", xlim = c(0,139), ylab = "Currently Hospitalized COVID Patients", main = "7 Day Currently Hospitalized Patients Forecast")
lines(seq(133,139,1), as.data.frame(currHospForecasts7)$fcst, type = "l", col = "red")

b. 90 Day

#90 day Forecasts
plot(seq(1,132,1), us$hospitalizedCurrently, type = "l", xlim = c(0,222), ylab = "Currently Hospitalized COVID Patients", main = "90 Day Currently Hospitalized Patients Forecast")
lines(seq(133,222,1), as.data.frame(currHospForecasts90)$fcst, type = "l", col = "red")

  1. Calculate ASE # WHERE YOU LEFT OFF
#varASE = mean((us$hospitalizedCurrently[123:132]-pred.var$fcst$y1[1:10])^2)
#varASE
  1. Windowed ASE

Multilayer Perceptron Model (MLP) / NN

Since we have been looking at additional regressors for all our above models and know that the best regressor to leverage will be positive increase of COVID cases. This regressor reflects similar behavior to that of current hospitlizations and can be model properly against hospitalizations. 1. Create Current Hospitalized ts variable

tUScurrhop = ts(us$hospitalizedCurrently)
  1. Create data frame of regressor: positive increase covid cases variable
tUSx = data.frame(positiveIncrease = ts(us$positiveIncrease))
  1. Forecast of positive increase of COVID cases
    1. 7 Day Forecast for Regressor
fit.mlp.new = mlp(ts(tUSx), reps = 50, comb = "mean")
mlp.new.fore7 = forecast(fit.mlp.new, h = 7)
mlp.new.fore7
##     Point Forecast
## 133       59916.13
## 134       63405.99
## 135       65059.79
## 136       66164.89
## 137       65396.11
## 138       64985.44
## 139       64589.15
b. 90 Day Forecast for Regressor
mlp.new.fore90 = forecast(fit.mlp.new, h = 90)
mlp.new.fore90
##     Point Forecast
## 133       59916.13
## 134       63405.99
## 135       65059.79
## 136       66164.89
## 137       65396.11
## 138       64985.44
## 139       64589.15
## 140       64839.67
## 141       64978.63
## 142       65158.69
## 143       65112.47
## 144       65094.63
## 145       65037.61
## 146       65056.91
## 147       65061.11
## 148       65084.50
## 149       65079.34
## 150       65081.11
## 151       65073.65
## 152       65077.61
## 153       65078.03
## 154       65082.25
## 155       65080.88
## 156       65080.74
## 157       65078.71
## 158       65079.45
## 159       65079.87
## 160       65081.28
## 161       65081.23
## 162       65081.13
## 163       65080.39
## 164       65080.36
## 165       65080.41
## 166       65080.86
## 167       65080.96
## 168       65080.99
## 169       65080.78
## 170       65080.73
## 171       65080.70
## 172       65080.81
## 173       65080.86
## 174       65080.88
## 175       65080.83
## 176       65080.81
## 177       65080.80
## 178       65080.83
## 179       65080.84
## 180       65080.85
## 181       65080.84
## 182       65080.83
## 183       65080.83
## 184       65080.83
## 185       65080.84
## 186       65080.84
## 187       65080.84
## 188       65080.84
## 189       65080.83
## 190       65080.84
## 191       65080.84
## 192       65080.84
## 193       65080.84
## 194       65080.84
## 195       65080.84
## 196       65080.84
## 197       65080.84
## 198       65080.84
## 199       65080.84
## 200       65080.84
## 201       65080.84
## 202       65080.84
## 203       65080.84
## 204       65080.84
## 205       65080.84
## 206       65080.84
## 207       65080.84
## 208       65080.84
## 209       65080.84
## 210       65080.84
## 211       65080.84
## 212       65080.84
## 213       65080.84
## 214       65080.84
## 215       65080.84
## 216       65080.84
## 217       65080.84
## 218       65080.84
## 219       65080.84
## 220       65080.84
## 221       65080.84
## 222       65080.84
  1. Combine observed new cases + forecast new cases
    1. 7 Day regressor var
new.regressor7 <- data.frame(c(us$positiveIncrease, mlp.new.fore7$mean))
new.regressor7 
##     c.us.positiveIncrease..mlp.new.fore7.mean.
## 1                                      3613.00
## 2                                      3171.00
## 3                                      4671.00
## 4                                      6255.00
## 5                                      6885.00
## 6                                      9259.00
## 7                                     11442.00
## 8                                     10632.00
## 9                                     12873.00
## 10                                    17656.00
## 11                                    19044.00
## 12                                    19724.00
## 13                                    19655.00
## 14                                    21897.00
## 15                                    24694.00
## 16                                    25773.00
## 17                                    27992.00
## 18                                    31945.00
## 19                                    33252.00
## 20                                    25550.00
## 21                                    28937.00
## 22                                    30737.00
## 23                                    30534.00
## 24                                    34419.00
## 25                                    34351.00
## 26                                    30561.00
## 27                                    27844.00
## 28                                    25133.00
## 29                                    25631.00
## 30                                    30298.00
## 31                                    30923.00
## 32                                    31964.00
## 33                                    27963.00
## 34                                    27468.00
## 35                                    25872.00
## 36                                    26333.00
## 37                                    28916.00
## 38                                    31784.00
## 39                                    34174.00
## 40                                    35901.00
## 41                                    27380.00
## 42                                    22605.00
## 43                                    25219.00
## 44                                    26641.00
## 45                                    29568.00
## 46                                    33056.00
## 47                                    29185.00
## 48                                    25767.00
## 49                                    22523.00
## 50                                    22449.00
## 51                                    25056.00
## 52                                    27490.00
## 53                                    27605.00
## 54                                    24810.00
## 55                                    21504.00
## 56                                    18236.00
## 57                                    22663.00
## 58                                    21193.00
## 59                                    26705.00
## 60                                    24710.00
## 61                                    24701.00
## 62                                    20171.00
## 63                                    20992.00
## 64                                    20853.00
## 65                                    21319.00
## 66                                    26513.00
## 67                                    24555.00
## 68                                    21590.00
## 69                                    20105.00
## 70                                    18798.00
## 71                                    16629.00
## 72                                    19385.00
## 73                                    22601.00
## 74                                    23522.00
## 75                                    23762.00
## 76                                    21639.00
## 77                                    20415.00
## 78                                    19982.00
## 79                                    20312.00
## 80                                    20839.00
## 81                                    23352.00
## 82                                    23106.00
## 83                                    18823.00
## 84                                    17012.00
## 85                                    17175.00
## 86                                    20803.00
## 87                                    22032.00
## 88                                    23477.00
## 89                                    25341.00
## 90                                    21373.00
## 91                                    18510.00
## 92                                    23435.00
## 93                                    23873.00
## 94                                    27527.00
## 95                                    31046.00
## 96                                    31994.00
## 97                                    27284.00
## 98                                    27017.00
## 99                                    33021.00
## 100                                   38684.00
## 101                                   39072.00
## 102                                   44421.00
## 103                                   43783.00
## 104                                   41857.00
## 105                                   39175.00
## 106                                   47462.00
## 107                                   50674.00
## 108                                   53826.00
## 109                                   54223.00
## 110                                   54734.00
## 111                                   45789.00
## 112                                   41600.00
## 113                                   51766.00
## 114                                   62147.00
## 115                                   58836.00
## 116                                   66645.00
## 117                                   63007.00
## 118                                   60978.00
## 119                                   58465.00
## 120                                   62879.00
## 121                                   65382.00
## 122                                   70953.00
## 123                                   77233.00
## 124                                   65180.00
## 125                                   64884.00
## 126                                   56971.00
## 127                                   63642.00
## 128                                   69150.00
## 129                                   71027.00
## 130                                   75193.00
## 131                                   65413.00
## 132                                   61713.00
## 133                                   59916.13
## 134                                   63405.99
## 135                                   65059.79
## 136                                   66164.89
## 137                                   65396.11
## 138                                   64985.44
## 139                                   64589.15
b. 90 day regressor var
new.regressor90 <- data.frame(c(us$positiveIncrease, mlp.new.fore90$mean))
new.regressor90
##     c.us.positiveIncrease..mlp.new.fore90.mean.
## 1                                       3613.00
## 2                                       3171.00
## 3                                       4671.00
## 4                                       6255.00
## 5                                       6885.00
## 6                                       9259.00
## 7                                      11442.00
## 8                                      10632.00
## 9                                      12873.00
## 10                                     17656.00
## 11                                     19044.00
## 12                                     19724.00
## 13                                     19655.00
## 14                                     21897.00
## 15                                     24694.00
## 16                                     25773.00
## 17                                     27992.00
## 18                                     31945.00
## 19                                     33252.00
## 20                                     25550.00
## 21                                     28937.00
## 22                                     30737.00
## 23                                     30534.00
## 24                                     34419.00
## 25                                     34351.00
## 26                                     30561.00
## 27                                     27844.00
## 28                                     25133.00
## 29                                     25631.00
## 30                                     30298.00
## 31                                     30923.00
## 32                                     31964.00
## 33                                     27963.00
## 34                                     27468.00
## 35                                     25872.00
## 36                                     26333.00
## 37                                     28916.00
## 38                                     31784.00
## 39                                     34174.00
## 40                                     35901.00
## 41                                     27380.00
## 42                                     22605.00
## 43                                     25219.00
## 44                                     26641.00
## 45                                     29568.00
## 46                                     33056.00
## 47                                     29185.00
## 48                                     25767.00
## 49                                     22523.00
## 50                                     22449.00
## 51                                     25056.00
## 52                                     27490.00
## 53                                     27605.00
## 54                                     24810.00
## 55                                     21504.00
## 56                                     18236.00
## 57                                     22663.00
## 58                                     21193.00
## 59                                     26705.00
## 60                                     24710.00
## 61                                     24701.00
## 62                                     20171.00
## 63                                     20992.00
## 64                                     20853.00
## 65                                     21319.00
## 66                                     26513.00
## 67                                     24555.00
## 68                                     21590.00
## 69                                     20105.00
## 70                                     18798.00
## 71                                     16629.00
## 72                                     19385.00
## 73                                     22601.00
## 74                                     23522.00
## 75                                     23762.00
## 76                                     21639.00
## 77                                     20415.00
## 78                                     19982.00
## 79                                     20312.00
## 80                                     20839.00
## 81                                     23352.00
## 82                                     23106.00
## 83                                     18823.00
## 84                                     17012.00
## 85                                     17175.00
## 86                                     20803.00
## 87                                     22032.00
## 88                                     23477.00
## 89                                     25341.00
## 90                                     21373.00
## 91                                     18510.00
## 92                                     23435.00
## 93                                     23873.00
## 94                                     27527.00
## 95                                     31046.00
## 96                                     31994.00
## 97                                     27284.00
## 98                                     27017.00
## 99                                     33021.00
## 100                                    38684.00
## 101                                    39072.00
## 102                                    44421.00
## 103                                    43783.00
## 104                                    41857.00
## 105                                    39175.00
## 106                                    47462.00
## 107                                    50674.00
## 108                                    53826.00
## 109                                    54223.00
## 110                                    54734.00
## 111                                    45789.00
## 112                                    41600.00
## 113                                    51766.00
## 114                                    62147.00
## 115                                    58836.00
## 116                                    66645.00
## 117                                    63007.00
## 118                                    60978.00
## 119                                    58465.00
## 120                                    62879.00
## 121                                    65382.00
## 122                                    70953.00
## 123                                    77233.00
## 124                                    65180.00
## 125                                    64884.00
## 126                                    56971.00
## 127                                    63642.00
## 128                                    69150.00
## 129                                    71027.00
## 130                                    75193.00
## 131                                    65413.00
## 132                                    61713.00
## 133                                    59916.13
## 134                                    63405.99
## 135                                    65059.79
## 136                                    66164.89
## 137                                    65396.11
## 138                                    64985.44
## 139                                    64589.15
## 140                                    64839.67
## 141                                    64978.63
## 142                                    65158.69
## 143                                    65112.47
## 144                                    65094.63
## 145                                    65037.61
## 146                                    65056.91
## 147                                    65061.11
## 148                                    65084.50
## 149                                    65079.34
## 150                                    65081.11
## 151                                    65073.65
## 152                                    65077.61
## 153                                    65078.03
## 154                                    65082.25
## 155                                    65080.88
## 156                                    65080.74
## 157                                    65078.71
## 158                                    65079.45
## 159                                    65079.87
## 160                                    65081.28
## 161                                    65081.23
## 162                                    65081.13
## 163                                    65080.39
## 164                                    65080.36
## 165                                    65080.41
## 166                                    65080.86
## 167                                    65080.96
## 168                                    65080.99
## 169                                    65080.78
## 170                                    65080.73
## 171                                    65080.70
## 172                                    65080.81
## 173                                    65080.86
## 174                                    65080.88
## 175                                    65080.83
## 176                                    65080.81
## 177                                    65080.80
## 178                                    65080.83
## 179                                    65080.84
## 180                                    65080.85
## 181                                    65080.84
## 182                                    65080.83
## 183                                    65080.83
## 184                                    65080.83
## 185                                    65080.84
## 186                                    65080.84
## 187                                    65080.84
## 188                                    65080.84
## 189                                    65080.83
## 190                                    65080.84
## 191                                    65080.84
## 192                                    65080.84
## 193                                    65080.84
## 194                                    65080.84
## 195                                    65080.84
## 196                                    65080.84
## 197                                    65080.84
## 198                                    65080.84
## 199                                    65080.84
## 200                                    65080.84
## 201                                    65080.84
## 202                                    65080.84
## 203                                    65080.84
## 204                                    65080.84
## 205                                    65080.84
## 206                                    65080.84
## 207                                    65080.84
## 208                                    65080.84
## 209                                    65080.84
## 210                                    65080.84
## 211                                    65080.84
## 212                                    65080.84
## 213                                    65080.84
## 214                                    65080.84
## 215                                    65080.84
## 216                                    65080.84
## 217                                    65080.84
## 218                                    65080.84
## 219                                    65080.84
## 220                                    65080.84
## 221                                    65080.84
## 222                                    65080.84
  1. Fit model for currently hospitalized w/ regressor of new cases
mlp.fit1 = mlp(tUScurrhop, xreg = tUSx, outplot = T, comb = "mean")

plot(mlp.fit1)

  1. Forecast w/ known Positive Increase Case Variable
    1. Currently Hospitalized 7 Day Forecast w/ Positive Increaser Regressor
currhosp.fore7 = forecast(mlp.fit1, h = 7, xreg = new.regressor7)
plot(currhosp.fore7)

b. Currently Hospitalized 90 Day Forecast w/ Positive Increaser Regressor

currhosp.fore90 = forecast(mlp.fit1, h = 90, xreg = new.regressor90)
plot(currhosp.fore90)

  1. Windowed ASE
  1. 7 Day Forecast
#if (!require(devtools)) {
#  install.packages("devtools")
#}
#devtools::install_github("josephsdavid/tswgewrapped", build_vignettes = TRUE)
library(tswgewrapped)
  • Train/ Test Data Sets
data_train.m <- data.frame(hospitalizedCurrently = us$hospitalizedCurrently[1:122], positiveIncrease = us$positiveIncrease[1:122])
data_test.m <- data.frame(hospitalizedCurrently = us$hospitalizedCurrently[123:132], positiveIncrease = us$positiveIncrease[123:132])
# search for best NN hyperparameters in given grid
model.m = tswgewrapped::ModelBuildNNforCaret$new(data = data_train.m, var_interest = "hospitalizedCurrently",
                                               search = 'random', tuneLength = 5, parallel = TRUE,
                                               batch_size = 50, h = 7, m = 7,
                                               verbose = 1)
## Registered S3 method overwritten by 'lava':
##   method     from   
##   print.pcor greybox
## Aggregating results
## Selecting tuning parameters
## Fitting reps = 10, hd = 2, allow.det.season = TRUE on full training set
## - Total Time for training: : 59.256 sec elapsed
  • The ASEs associated with the grid of hyperparameters is shown in the table and heatmap below.
res.m <- model.m$summarize_hyperparam_results()
res.m
##   reps hd allow.det.season     RMSE      ASE   RMSESD    ASESD
## 1   10  2             TRUE 3103.429 13970241 2184.689 16296534
## 2   11  3            FALSE 3129.341 14942705 2380.109 19673017
## 3   16  1            FALSE 3559.638 16480775 2047.127 15988977
## 4   16  3             TRUE 3204.752 15391189 2373.358 19509361
## 5   22  3            FALSE 3268.191 16196849 2463.200 19979406
model.m$plot_hyperparam_results()

  • Windowed ASE: The best hyperparemeters are shown below
best.m <- model.m$summarize_best_hyperparams()
best.m
##   reps hd allow.det.season
## 1   10  2             TRUE
final.ase.m <- dplyr::filter(res.m, reps == best.m$reps &
                    hd == best.m$hd &
                    allow.det.season == best.m$allow.det.season)[['ASE']]
final.ase.m
## [1] 13970241
  • The best hyperparameters based on this grid search are 13 repetitions and 1 hidden layers, and allow.det.season = FALSE which has a mean windowed ASE.
# Final Model
caret_model.m = model.m$get_final_models(subset = 'a')
caret_model.m$finalModel
## MLP fit with 2 hidden nodes and 10 repetitions.
## Series modelled in differences: D1.
## Univariate lags: (1,2,3)
## 1 regressor included.
## - Regressor 1 lags: (1,3)
## Forecast combined using the median operator.
## MSE: 1065689.5459.
  • Plot final model
# Plot Final Model
plot(caret_model.m$finalModel)